
#include "mex.h"
#include <string.h>
#include <math.h>
#include "cplx.h"
#include <stdio.h>

#define PMAX 20
#define CMAX 20

double erf(double y);
double rh1(double z, double P[], double C[]);
double beta1(double z, double P[], double C[]);
double step(double x);
double rect(double x, double d);
double calculateReflec(double P[], double C[], double x);



void mexFunction(int nlhs, mxArray *plhs[],
        int nrhs, const mxArray *prhs[]) {
    double P[PMAX+1], C[CMAX+1], *qz;
    double *p_ptr, *c_ptr;
    double *RF;
    int mrows, ncols;
    int n;
    
    p_ptr=mxGetPr(prhs[0]);
    P[0]=0;
    for (n=1;n<=PMAX;n++) {
        P[n]=p_ptr[n-1];
    }
    c_ptr=mxGetPr(prhs[1]);
    C[0]=0;
    for (n=1;n<=CMAX;n++) {
        C[n]=c_ptr[n-1];
    }
    
    qz=mxGetPr(prhs[2]);
    mrows=mxGetM(prhs[2]);
    ncols=mxGetN(prhs[2]);
    plhs[0]=mxCreateDoubleMatrix(mrows, ncols, mxREAL);
    RF=mxGetPr(plhs[0]);
    for ( n=0;n<mrows ;n++) {
        RF[n]=calculateReflec(P, C, qz[n]);
    }
}

/*This following reflectivity part is written by Wenjie Wang in March,2010*/
double calculateReflec(double P[], double C[], double x)       /* !!!--!!!*/
{
  /*fitref25_AA_APS.5.c
   Wenjie WANG Dec,2010
   */
  struct cplx f[203], an[203], rn[203], rn0, pj, numerator, denominator,temp ;/*pj is term in Parratt's recursion of layer N.
                                                                                      see book page 82*/;
  double rho[203],my[203], del[203],d[203], bet[203], y, d_slice;/*In this two-box case, the slice is only 2, so the array length is also not short*/
  double length;
  double slu, theta, frsnll,xnew,z; 

  int nmax, n, j, i;
  nmax      = C[1] + 1;
  rho[0]    = 0.0;
  my[0]     = 0.0;
  d[0]      =0;
  d[nmax]    =0;
  rho[nmax] = P[4];
  my[nmax]  = P[5];
 
 
  slu       = 1.0;  /* For neutrons. */
  
  if (C[7] == 0.0) slu = 2.817938e-5;     /* For X-rays. r-e in [AA] */

  if (C[2] == 0.0) {                            /* Refelectivity Calculations */
    xnew=x-P[1];
    theta = xnew/(2.0*C[5]);/* x is one single value here corresponding to qz */
    
    del[nmax] = slu * 2.0 * 3.14159 / (C[5]*C[5]) * rho[nmax];/* nmax ->This delta of subphase*/ 
    bet[nmax] = .5/C[5]*1.e-8*my[nmax];/* my is mu in [cm^-1]. bet is dimensionless. C[5] is Angstrom ^-1 */
    f[nmax] = cplxsqrt(makecplx(theta*theta - 2.0 * del[nmax],2.0 * bet[nmax]));/* Sign of Bet is Positive*/
    /* f[..] is moment transfer kz; sign in front of bet[] is contraversial */
    del[0] = 0.0;
    bet[0] = 0.0;
    
    /*  P[nnmax*4-4]=258/P[nnmax*4-6]/C[10];*/
    /*  P[nnmax*4]=0;*/
  
    f[0] = makecplx(theta,0); 
    
    /* Calculate Fresnel law for subphase-only: */
  
    rn0= cdiv( cdiff(f[0],f[nmax]),csum(f[0],f[nmax]) );
    frsnll=modulus(rn0)*modulus(rn0);
    y=frsnll;

    /* The following loop is asign the value to boxes
    Index 0-> air
    Index nmax->subphase
    Index nmax-1 -> the layer(box) above the subphase
    Index 1 -> the layer(box) bordered on air
     */

        
    z=-P[8]*5; 
    length=3*P[8]+P[6]+3*P[11]+P[9]+3*P[14]+P[12]+3*P[15]; /* Total Length */ 
    d_slice=length/C[1];
    
    for (n=nmax-1;n>=1;n--) {
      d[n]=d_slice;
    }

   
    for (n=nmax-1;n>=1;n--) {
     
      z=z+0.5*d[n+1]+0.5*d[n];
     
      rho[n]=rh1(z,P,C);
      my[n]=0;/*beta1(z,P,C);*//* The absorption is taken into account by call on beta()*/
      del[n]=slu*2.*3.14159/(C[5]*C[5])*rho[n];  /*Refraction component delta*/     
      bet[n]=.5/C[5]*1.e-8*my[n];                /*Refraction component beta */
      f[n]=cplxsqrt(makecplx(theta*theta-2.0*del[n],2.0*bet[n]));/* Sign of Bet is Positive*/
    }

    rn[nmax]=makecplx(0.0,0.0);/*Reflection from bottom of subphase is zero*/
       
    for (n=nmax; n>=1; n--) {
      
      an[n-1]=cdiv( cdiff(f[n-1],f[n]),csum(f[n-1],f[n]) );/*Ri in Weibu's dissertation eq.(2.44) */
      pj=cplxexp(cmult(makecplx(0,2.*C[5]*d[n]),f[n])); /* exp[2j*k(i+1)*d(i+1)]*/
      temp=cmult(rn[n],pj); /* term r(i+1)*pj **/
      numerator=csum(an[n-1],temp);/* R(j)+r(j+1)*pj */      
      denominator=csum(makecplx(1.0,0),cmult(an[n-1],temp));/*1+R(j)*r(j+1)*pj */
      rn[n-1]=cdiv(numerator,denominator);/*r[0] is what we are looking for*/
    }

/* Output Of Reflectivity: for unbound head group */
    
    y = P[3]*modulus(rn[0])*modulus(rn[0]);

    if (C[6] != 0.0) y = y/frsnll; 


    return y;
  } /* End of Reflectivity Calculations */

  /* Real Space (Rho vs. X=Height: */

  if (C[2] == 1) {   /*Electron Number Density Profile*/ 
    y=rh1(x,P,C);
    return y;
   }
}


double rh1(double z,double P[],double C[]) 
{
  double zi[4],rho[4],sig[4]; 
  double w0,w1,w2,w3,w4;                                /*weight*/
  double eta1,eta2,eta3;
  double rsum=0;

  
  zi[0]=0;                                 /*Subphase-Layer-1 interface*/
 
  zi[1]=P[6];                              /*Layer-1 Layer-2 interface*/

  zi[2]=zi[1]+P[9];                        /* Layer-2 Layer-3 interface*/
   
  zi[3]=zi[2]+P[12];                       /* Layer-3 Layer-4 interface */

  rho[0]=P[4];                             /*Electron density of subphase*/

  rho[1]=P[7];                             /*Electron density of Layer 1 */

  rho[2]=P[10];                             /*Electron density of Layer 2*/

  rho[3]=P[13];                             /* Electron density of Layer 3*/

  
  sig[0]=P[8];                             /* Roughness of subphase-Layer-1 interface */

  sig[1]=P[11];                            /* Roughness of Layer-1 & Layer-2 interface */
  
  sig[2]=P[14];                             /* Roughness of Layer-2 & Layer-3 interface */

  sig[3]=P[15];                            /* Roughness of Layer-3 & Layer-4 interface */

 
  eta1=(zi[0]*sig[1]+zi[1]*sig[0])/(sig[0]+sig[1]);

  eta2=(zi[1]*sig[2]+zi[2]*sig[1])/(sig[1]+sig[2]);

  eta3=(zi[2]*sig[3]+zi[3]*sig[2])/(sig[2]+sig[3]);

  w0=0.5*(1-erf((z-zi[0])/sig[0]));

  if (z <= eta1 ){
      w1=0.5*(1+erf((z-zi[0])/sig[0]));
     }
  else {
      w1=0.5*(1-erf((z-zi[1])/sig[1]));
     }
  
  if (z <= eta2 ){
      w2=0.5*(1+erf((z-zi[1])/sig[1]));
     }
  else {
      w2=0.5*(1-erf((z-zi[2])/sig[2]));
     }
     
  if (z <= eta3 ){
      w3=0.5*(1+erf((z-zi[2])/sig[2]));
     }
  else {
      w3=0.5*(1-erf((z-zi[3])/sig[3]));
     }

  w4=0.5*(1+erf((z-zi[3])/sig[3]));

  rsum=(w0*rho[0]+w1*rho[1]+w2*rho[2]+w3*rho[3])/(w0+w1+w2+w3+w4);
  return rsum;
}




double beta1(double z,double P[],double C[])/* beta is absorption for unbound head group*/
{
    return 0.0;
}


/****  THIS IS END OF CODE WRITTEN BY WENJIE WANG IN JAN.2010 ***/

/*****************************************************************************************
Error function definition from Abramowitz and Segun page 299 7.1.25 */


double erf(double y)
{
  double yx, tx, zx,sgn;
  if (y>=0) sgn=1.0;
  else sgn=-1.0;

  yx=sqrt(y*y)/sqrt(2.);
  tx=1./(1. + .47047*yx);
  zx= 1.-(.34802*tx-.0958798*tx*tx+.7478556*tx*tx*tx)*exp(-yx*yx); 
  zx=sgn*zx;
  return zx;      
}
double step(double x)
{
  if (x>=0) return 1.0;
  else return 0.0;
}

double rect(double x,double d)
{
  if (x>=-d/2.0 && x< d/2.0) return 1.0;
  else return 0.0;
}
